Intrinsic noise induced resonance in presence of sub-threshold signal in Brusselator 
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In a system of non-linear chemical reactions called the Brusselator, we show that intrinsic noise 
can be regulated to drive it to exhibit resonance in the presence of a sub-threshold signal. The 
phenomena of periodic stochastic resonance and aperiodic stochastic resonance, hitherto studied 
mostly with extrinsic noise, is demonstrated here to occur with inherent systemic noise using exact 
stochastic simulation algorithm due to Gillespie. The role of intrinsic noise in a couple of other 
phenomena is also discussed. 
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Studies of noise driven regularity in non-linear 
systems have not been done as extensively for 
"intrinsic noise" as for extrinsic noise. Here we 
give direct demonstration of stochastic resonance 
(both periodic and aperiodic) in a chemical sys- 
tem with respect to intrinsic systemic fluctua- 
tions, using exact stochastic simulation method 
due to Gillespie. Moreover, an interplay of the 
intrinsic and extrinsic noises is analyzed for these 
noise invoked resonances. 



I. INTRODUCTION 

Chemical and biological systems quite generally are 
stochastic in nature. The stochasticity may arise from 
the inherent probabilistic nature of the processes involved 
or external environmental interferences — accordingly 
they are referred to as "intrinsic" and "extrinsic" noise, 
respectively. For a chemical system specified by certain 
reaction rate constants, the chemical reactions in actual 
reality vary randomly, and molecular numbers do fluc- 
tuate — this is intrinsic to the system tlL 2 [ . Similarly 
biological processes like gene expression [3|, |4| involve in- 
trinsic stochasticity. As opposed to this, if fluctuations 
in the medium or other chemical components external to 
the subsystem of interest, indirectly affect the biochem- 
ical processes involved, they are regarded as extrinsic. 
Although for biochemical systems with large number of 
molecules, deterministic approximations of chemical re- 
actions are often enough, for systems with small num- 
ber of molecules, intrinsic fluctuations are relatively large 
and stochastic treatment is necessary. In the latter case, 
the system's response may go beyond mere small excur- 
sions about the averages and unexpected behavior may 
appear due to random crossing of thresholds. A theoret- 
ical framework within which intrinsic noise in biochem- 
ical systems is studied is the chemical master equation 
(CME) approach. A widely used numerical method to 
exactly implement the CME is Gillespie algorithm [5- 7] , 
as the CME is often not easily analytically tractable. 

Noise, intrinsic or extrinsic, is commonly thought of 
as an undesirable disturbance. Yet in the realm of non- 
linear systems, it is well known now that noise can aid 
the output attain regularity, or at times, reveal hidden 
order. Here by "regularity" we imply either enhanced pe- 



riodicity of the output response, or enhanced correlation 
of the output response to the input signal. In this paper, 
we study primarily periodic stochastic resonance (PSR) 
and aperiodic stochastic resonance (ASR). These effects 
have been seen in multiple non-linear systems [8l-[l6| with 
extrinsic noise. When a non-linear system is near a bi- 
furcation threshold, and subjected to an weak input sig- 
nal superposed with a noise, the response of the output 
shows maximum "regularity" for an optimum strength of 
noise. The input signal used for PSR is periodic and for 
ASR is aperiodic. The superimposed noise used in earlier 
studies were mostly extrinsic — for example, in the the- 
oretical models [§-[lCjj noise has been added externally to 
the dynamical equations, while in experiments on elec- 
trochemical cell [H, EH external noise has been added to 
the voltage. 

A natural question is that can these phenomena be 
seen as a result of variation of intrinsic noise inherent 
to the system, near a bifurcation threshold? Experimen- 
talists vary different sources of internal noise to study 
regularity of response, suitable and specific to their sys- 
tem of interest — e.g. in an electrochemical cell chloride 
ion concentration was varied 

[13, 

while in hippocampal 
CA1 neurons sub-threshold cathodic current was used 
[H ] . Although sources of internal noise can be diverse in 
practical systems, for keeping our theoretical discussion 
general, we follow the standard ideas of intrinsic stochas- 
ticity in chemical reactions used by Gillespie d, Q . The 
idea is to simulate CME numerically using exact stochas- 
tic simulation algorithm given by Gillespie p| H|. In- 
stead of the latter direct approach, in the literature, often 
chemical Langevin equation (CLE) @, [H| has been used 
[20M22I I22I [23| | . It is important to note that CLE is an 
approximation of CME [24| , and sometimes produce mis- 
leading results [25[ . Therefore, CME is preferable. In this 
paper we implement CME using the Gillespie method. 

The phenomenon of coherence resonance (CR) p6l - [3C| , 
where there is no input sub-threshold signal, has been 
studied earlier using Gillespie and CLE method [20l - [22l 
l3ll - [34j |. To our knowledge, demonstration of PSR and 
ASR in a chemical system with Hopf bifurcation using 
exact CME approach is not well known; although this has 
been studied for bistable nonlinear chemical systems (35l . 
|36[. In the present work, we study PSR and ASR in the 
chemical Brusselator system using Gillespie algorithm. 
Furthermore, we explore the effect on PSR and ASR due 
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to the interplay of extrinsic and intrinsic noises. Apart 
from that, we revisit the CR and another noise driven 
phenomenon to make some interesting observations. 

The Brusselator oscillator 0, [23|, |37j system involves 
the following four chemical reactions. 

Z x ^ X :Zx fixed 
Z 2 + X H Y + Z 3 : Z 2 fixed 
2X + Y H 3X 



X % Zi 



(1) 



Zx and Z 2 are kept fixed by having infinite sources of 
these reactants. In Eq. [TJ {kj} (j — 1,2,3,4) are the de- 
terministic reaction rate constants for the four reactions 
respectively. The symbols X, Y. Zx, Z 2 , Z 3 , Zx repre- 
sent the molecular numbers for the six chemical species 
involved. The system volume is denoted by fi. The num- 
bers of X, Y, Z 3 and Zx vary with time; out of them X 
and Y are the primary variables of interest. For certain 
choice of the reaction rate parameters, the reactants X 
and Y fluctuate around fixed point values, while other 
choice of rates, they fluctuate around limit cycle behav- 
ior. These two regimes are separated by a Hopf bifur- 
cation. Within Gillespie algorithm, propensity functions 
@, @ for the above reactions are, 



a-3 



fii = Zxkx, a 2 = Z 2 Xk 2 /fl, 
X{X - l)Yk 3 /n 2 , and 5 4 = Xk A . 



(2) 



respectively. 

We discuss the Gillespie algorithm 0-0] for this sys- 
tem briefly. As Z\ and Z 2 , coming from infinite sources, 
are not time varying, we consider a state vector s(t) ex- 
cluding them, namely s(t) = (X(t), Y(t), Z 3 (t), Z A (t)). 
Starting with a state vector s(t) at time t, the j th reaction 
out of the four in Eq. [5]is chosen with a probability Oj/eto, 
after a waiting time St, drawn from a probability distri- 
bution function P(5t\s(t)) = a exp(— aoSt). Here {aj} is 
given by Eq. [3J and a = Ylj=x ano ^ these are func- 
tions of s(i). After the implementing j th reaction, the 
state vector s(t) is updated to s(t + St) =s(t)+Uj, where 
using Eq. [T]one can easily see that the state change vec- 
tors, v x = (1,0,0,0), v 2 = (-1,1,1,0), v 3 = (1,-1,0,0), 
and v± = (-1,0,0,1). 

The way to tune intrinsic noise strength to study reso- 
nances for this system, is to vary fi, Zx and Z 2 , keeping 
the concentrations zx = Zi/fl and z 2 = Z 2 /Q fixed by 
hand. Under the latter condition, we observed that the 
numbers X, Y also spontaneously adjust to keep their 
concentrations x = X/Q and y = Y/VL, on an average 
close to a constant. The observation is expected since 
the concentration x and y, if approximated to be deter- 
ministic, satisfy the following equations: 



dx 
~dt 
(kj_ 
dt 



k\Z\ — k 2 z 2 x + k 3 x 2 y — k^x 

k 2 z 2 x — k 3 x 2 y. (3) 
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FIG. 1: Shows the limit cycle regime (LCR) and the fixed 
point regime (FPR) in 'a-b' parameter space for Brusselator. 
Dynamical behavior of the points P, Q, R, and, S in the de- 
terministic and the stochastic cases are shown in Fig. [S] 



The above equations show that if the concentration z\ 
and z 2 are fixed, the concentration x and y can be deter- 
mined. Thus by increasing (decreasing) i7, Z\ and Z 2 , 
keeping concentrations zx and z 2 fixed, we ensure that we 
stay at the same parameter set point of the non-linear 
system, and keep decreasing (increasing) the "intrinsic 
noise" around it. The condition for getting limit cycle 
behavior for the deterministic Eq. [3] is (23j . 



b < (2a- 1)(1 - a) 



(4) 



where a = k 2 z 2 j (k 2 z 2 + k^) and b — (k 2 z 2 k 3 )/[(k 2 z 2 + 
kx) 3 )). In Fig. [U the solid line is Hopf bifurcation line 
separating fixed point regime (FPR) from the limit cycle 
regime (LCR). 

Although our actual simulation is stochastic following 
Eq. [l] in order to identify the average position of the set 
point of the system, we fall back whenever necessary to 
the Eq. [3] By this we mean that, a particular choice of 
kx, k 2 , k 3 , kx, Zx and z 2 , corresponds to a particular a 
and b value, and the Fig. [T] tells us the average location 
of the set point of the system. Thus using Fig. [T] as 
guidance, we proceed to check the phenomena of PSR 
and ASR. 



II. RESULTS 



In this section we show that the phenomena PSR and 
ASR happen due to intrinsic noise of the Brusselator sys- 
tem. In each case our set point is in the FPR, represented 
schematically by S shown in Fig. [T] We tune the intrin- 
sic noise strength for these reactions by varying Q, Z\ 
and Z 2l such that the concentrations zx = Zx/Sl and 
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FIG. 2: Data for PSR. (a): The lower curves in each box show behavior of x(t) for increasing system volumes (decreasing noise) 
: Q, — 0.02 (I), Q, — 0.14 (II), and £7 = 1.0 (III), in the presence of sub-threshold periodic signal s p (t) (the upper curves in every 
box). Here the pulse amplitude is -0.010, the interpulse interval is 10.0, and pulse width is 2.0. We observe that output spikes 
become more correlated with the input periodic signal in case II, compared to I and III. Note that an output spike begins with 
a reduction from the steady value, followed by a rapid rise — thus there is no significant phase shift between an input pulse 
and the output spike, as may appear at a first glance. Panel (b): The solid line with points shows absolute value of power 
norm Co against f2 (pure intrinsic noise); the dashed line shows the same, for an additional extrinsic noise (see the text). Note 
that | Co | has a peak at a smaller value of Q for the added extrinsic noise. 



Z2 = Z2/O, are held fixed by hand as discussed above. 
The noise makes the system occasionally excurse into 
LCR from FPR and one sees spikes in x and y as a re- 
sult of that. For different noise strengths we have looked 
at the output response x and y, but in the subsequent 
discussion we will focus only on the variable x for brevity. 

To keep the discussion simple we will suppress the 
explicit dimensions of the rate constants {kj}, and the 
volume n. This makes us avoid putting any units for 
x(t), t, Vn , and Co in all the subsequent figures 2- 
5. Furthermore, we denote the strength of the intrinsic 
noise by ft — the larger the £1, the smaller the intrinsic 
noise. Thus in all the figures 2-4, noise decreases along 
the increasing O axis. 

To characterize the regularity in the response of the 
system, a measure of the correlation between output re- 
sponse and input signal, called the power norm Cq fl3l . [lr| 
is used: 

Co = ((x(t) (x))(s(t) - («))). (5) 

Here (.) denotes the time average. The variable x(t) is 
the time varying output of the system at time t, and s(t) 
is the input signal. The signal s(t) is periodic for PSR 
and aperiodic for ASR. Co can be positive or negative 
depending on the relative sign of output response and 
the input signal. In this study Co is always positive, but 



for easier comparison with earlier publications [13l . Il6| we 
use |C |. 

We will decide the average location of the set point 
of our system, in 'a-b' space by choosing appropriate 
ki, k2, fe, ki, Zi, and Z2 values. In our simulation 
of PSR and ASR, discussed below, we take &2 = 0.1, 
k 3 = 0.00005, z x = 100 and z 2 = 500. The values of k x 
and &4 are different for different cases and are mentioned 
in the respective subsections. 

A. Periodic Stochastic resonance 

To study the PSR, we choose k\ and £4 to be time 
dependent: k\{t) — k 2 Z2(l — — s p(t))/( a o + s p(t)) an d 
/c 4 (i) = [b /(z1k 3 (ao + Sp(t)))] 1 / 2 . Here s p (t) is a time 
varying periodic signal. By doing periodic modulation of 
fei and &4 with time, the average location (a, b) of the 
set point gets modulated periodically with time and is 
given as (arj + s p (t),bo) — see the relations just below 
Eq. |U Thus s p (t) is to be regarded as a periodic input 
to the system. We take arj = 0.980, bo = 0.001, and 
the amplitude of s p is -0.010, such that minimum value 
of a = 0.980 - 0.010 = 0.970, is still greater than the 
threshold value a t h = 0.967 of the Hopf bifurcation point 
{dth, bo) in the absence of noise. The sub-threshold signal 
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FIG. 3: Data for ASR. Panel (a): The lower curves in each box show behavior of x(t) for increasing system volumes (decreasing 
noise) : 17 = 0.02 (I), 17 = 0.15 (II), and 17 = 1.40 (III), in the presence of sub-threshold aperiodic signal s a (t) (the upper 
curves in every box). Here the pulse amplitude is -0.010, the pulse width is 2.0. Interpulse intervals are aperiodic and its exact 
description is given in the text. We observe that output spikes become more correlated with the input aperiodic signal in case 
II, compared to I and III. As commented in Fig. [2l there is no significant phase shift between an input pulse and the output 
spike, as may appear at a first glance. Panel (b): The solid line with points shows absolute value of power norm Co against 17 
(pure intrinsic noise); the dashed line shows the same, for an additional extrinsic noise (see the text). Note that in this case 
also | Co | has a peak at a smaller value of 17 for the added extrinsic noise. 



amplitude value is 77 percentage of the threshold value 
(separation between the average set point position and 
bifurcation point). For these parameters and the input 
signal we run our Gillespie simulation for different val- 
ues of 17, Z±, and Z%, thereby varying the intrinsic noise 
strength. Fig. [5(a) shows output response (x variable) 
for three 17 values — (I) low, (II) intermediate and (III) 
high. We see that for an intermediate 17 (II), the output 
spike pattern becomes more correlated with the input pe- 
riodic pulse pattern, as compared to the cases I and III 
- this is the point of resonance. From the time series 
x(t) and s p {t), we calculate Cq using Eq. [5] The solid 
line with points in Fig. [2(b) shows \Cq\ against 17, in the 
absence of any extrinsic noise. An unimodal behavior 
is seen, where the maximum (implying maximal regu- 
larity) corresponds to the optimum value of 17 = 0.14. 
Thus with intrinsic noise we have demonstrated PSR in 
the Brusselator system. 

If extrinsic noise was present in the system, would the 
internal noise causing the phenomenon of PSR, act to 
reinforce or subdue the former effect? To explore this 
interesting question of interplay between extrinsic and 
intrinsic noises, we proceed to add a noise term to the 
average set point position: a — > a + D£. We choose 
extrinsic noise £ to be a uniform box distribution between 
[0, 1] and D = 0.002. Thus the noise is one sided and 



pushes (randomly though) the average set point further 
away from the bifurcation point. The question is whether 
we will require larger or smaller internal noise to achieve 
PSR in the presence of finite D. The dashed curve in Fig. 
[2(b) shows that the resonance occurs at a smaller value 
of 17 (i.e. larger intrinsic noise). Thus we observe that 
the two sources of noise, extrinsic and intrinsic, actually 
superimpose. 



B. Aperiodic Stochastic resonance 

ASR is similar to PSR, except that the signal s a is 
aperiodic. The input aperiodic signal used by us is an 
aperiodic pulse train. Aperiodicity comes only in the in- 
terpulse interval of the signal, which given by 6.0+8.0 x r. 
Here r is a random number drawn from uniform box dis- 
tribution between to 1. Note that the width of the 
pulses are the same as s p . We take the same parameters 
as in PSR, except the time varying k±(t) = fc2Z2(l — &a — 
s a (t))/(a + s a (t)) and h(t) = [b /(zfk 3 (a + s a (t)))}^ 2 . 
The pulse amplitude remains the same as for s p (t) (see 
section III A[) such that the average set point location re- 
mains greater than a t h = 0.967. Again we vary the in- 
trinsic noise by varying 17, Z\, and Z2. The Fig. [3(a) 
shows output response [x variable) for three values of 17: 
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(I) high, (II) intermediate and (III) large. We see that 
for an intermediate Q (II) output spike pattern become 
more correlated with the input aperiodic pulse pattern as 
compared to I and III — this is the point of resonance. 
From the time series we calculate Co using Eq. [5j and the 
solid line with points in Fig. EJb) shows |Co| against ft, 
in the absence of extrinsic noise. The maximum regular- 
ity corresponds to the optimum fl = 0.15, just like PSR. 
Thus with intrinsic noise we have demonstrated ASR in 
the Brusselator system. 

Analogous to PSR, in this case too we explore the in- 
terplay between extrinsic and intrinsic noises. We add 
a noise term to the average set point position: ag — > 
ao + D£. Extrinsic noise £ is uniformly distributed over 
[0, 1] and D = 0.002. The dashed curve in Fig. ^b) 
shows that the resonance occurs at a smaller value of Vl 
(i.e. larger intrinsic noise). Thus for ASR, we also ob- 
serve that the extrinsic and intrinsic noises superimpose. 
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FIG. 4: Data for CR. Normalized variance Vn plotted against 
fi — data from calculation using time series and Eq. [6] (sym- 
bols A joined by line) alongside data obtained using Eq. [7] 
(a). Vn has a minimum for an optimum value of Q, = 0.25. 



C. Test of a mathematical formula for coherence 
resonance 

As noted, coherence resonance (CR) has been studied 
earlier with intrinsic noise using exact stochastic simu- 
lations [3ll - [34j |. In this sub-section, we revisit this phe- 
nomenon to test a semi-analytical formula for normalised 
variance Vn recently proposed [3(| ■ In the case of CR 
there is no input signal, i.e s(t) = 0, and so Co (Eq. 
[5]) cannot be used as a measure of regularity. Instead 
normalized variance Vjy [HI, H3| is used: 

Vn = VVp 2 > - (r P f/(r P ), (6) 

where t p is interspike time intervals and (.) again denotes 
time average. Typically Vn is enumerated from the time- 
series analysis of spikes generated by the system, using 
Eq. [BJ An alternative formula for Vn as a function of 
noise strength was proposed recently, which directly re- 
flects the theoretical understanding of the phenomenon 
of CR. The formula 30] is as follows: 

V = T CSc(^) ,„s 

N r min (n) + T esc (n)' {,) 

Unlike Eq. O this is not a definition of Vn , but a sug- 
gested mathematical form. It expresses the important 
fact that CR comes about due to competitive interplay 
of two time scale t csc (a first passage activation time) 
and T m i n (an excursion time of virtual limit cycle) [30| . 
The formula was tested in [3(| for the FitzHung-Nagumo 
system (2(| and a chemical oscillator model [38[ in the 
presence of extrinsic noise. Here we would like to see if it 
holds good for intrinsic noise. Note again that intrinsic 
noise strength is proportional to £2 , the analog of ex- 
trinsic noise strength D [30]. In Eq. [7]t csc and r min can 
be calculated from the probability distribution function 
P(r p ) of T p . The function P{t p ) remains zero between 



T p = and t p — r ro j n , and then rises to a peak and even- 
tually decays with an exponential tail for large t p with a 
time constant T esc . 



We choose k\ = 16.551 and k^ = 1.546 which corre- 
spond to a = 0.970 and b = 0.001, setting the average 
location of the set point of our system in the FPR. Vn 
was calculated from the time series analysis of the output 
signal x(t) for various Vl — the curve is plotted in Fig. [4] 
(with A symbols joined by a continuous line). Next, we 
calculated P{t p ) and then found r osc and r m ; n from it — 
we used the latter to calculate Vn using Eq. [7] (shown 
with ▲ symbols in Fig. 0]). The curves of Vn obtained 
by two different procedures match reasonably well con- 
firming the validity of Eq. in the case of intrinsic noise 
in Brusselator. 



D. Output response on varying reactant molecular 
numbers, keeping Q fixed. 

In certain situations, it may be inconvenient to keep 
reactant concentrations fixed. It may be convenient to 
vary just the reactant numbers, or just the volume, in- 
dependently. Does any interesting phenomenon arise in 
such cases? Here we consider varying numbers Z\ and 
Zi, keeping Q fixed. This makes average concentrations 
of the reactants Z\, Z2, x, and y vary, and hence make the 
average set point of the system drift. The drift carries the 
system from FPR through LCR, and the noisy observable 
x(t) shows interesting change in temporal pattern due to 
that. The latter evolution is explicitly shown in Fig. [Ha). 
By tuning Z\ from 888 through 19126, and Z-z from 11 
through 419 (with fl = 1.0), the x(t) behavior changes 
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FIG. 5: Shows behavior of x(t) in four different points P, Q, R, and S (from bottom to top) of 'a-b' space. Panel (a) are the 
output responses in the presence of intrinsic noise and panel (b) are the corresponding deterministic responses. For these data, 
we used ki = 0.10, k 2 = 0.10, fc 3 = 0.00005, fc 4 = 1.75, and Q = 1.0. The various values of (Zi, Z 2 ) pairs used are (888, 11) 
for P, (1270, 21) for Q, (8421, 201) for R, (19126, 419) for S. 



from the bottom to the top of Fig. [SJa). The spiking 
behavior of x(t) passes successively from irregular— > reg- 
ular — > irregular, reminding us of the phenomenon of CR. 
But it should be noted that this phenomenon has noth- 
ing to do with CR. The behavior seen in Fig. [5ja) are 
noisy excursions on top of the steady behavior at P, Q, 
R, and, S (see Fig. [1]), respectively. Due to the change 
of concentrations, the system drift from FPR— !> LCR— !> 
FPR. To substantiate the latter assertion, we show the 
deterministic steady behavior at P, Q, R, and S in Fig. 
EJb) corresponding to their counterparts in Fig. EJa) - 
the similarity is obvious. Thus the above way of varying 
reactant concentrations keeping D, fixed, not only varies 
internal noise, but also make the set point of the system 
drift along with it. The effect is interesting but this pro- 
cedure is unsuitable for tuning inherent noise to study 
resonance phenomena in the system. 



system can induce enhanced "regularity" in response, 
seen in phenomena like PSR and ASR, just like extrin- 
sic noise. The intrinsic noise strength has been varied 
for given set {%}, by varying system volume fi and re- 
actants numbers, keeping the concentration of the reac- 
tants z\ and Z2 fixed; the average concentrations of x 
and y spontaneously adjusted to stay constant. In sec- 
tion III D[ we discussed that if we vary the number of 
the reactants keeping f2 fixed, their concentrations vary, 
leading to drift of the average set point of the system. 
Interesting variation of temporal pattern may follow due 
to drift from FPR— s> LCR— s> FPR, but that has nothing 
to do with the phenomenon of CR. In section Hi CI a sim- 
ple formula for the normalised variance Vjv as function 
of Q, (inverse strength of intrinsic noise) was tested. 



III. DISCUSSION 



In this paper, we have shown that an optimal value 
of the intrinsic noise present in a chemical Brusselator 
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